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This paper studies charge distribution in buried-channel charge- 
coupled devices. Detailed development of a one-dimensional electrostatic 
model is presented and a numerical solution of the resulting nonlinear po- 
tential equations is described. Graphical results show the charge-filling 
mechanism and the relationship between the oxide-semiconductor interface 
potential and total free positive charge. 

I. INTRODUCTION 

This paper describes a numerical determination of the distribution 
of charge in a one-dimensional model of a buried-channel charge- 
coupled device (CCD). 1 Several calculations have recently been made 
of the static potential in CCD's in the absence of stored charge. 2-4 In 
addition, the motion of the stored charge under dynamic conditions 
has been studied by means of essentially one-dimensional models which 
do not involve a true knowledge of the distribution of stored charge or 
the charge-carrying capacity of the CCD. 2 - 5 However, so far it has not 
been possible to calculate even the static stored-charge distribution in 
a two-dimensional model of a buried-channel CCD, much less to follow 
the motion of this charge under dynamic conditions. In this paper a 
start is made on the problem by calculating the static distribution of 
stored charge in a one-dimensional model of a buried-channel CCD. 
The resulting information on the charge distribution is of interest in 
itself. However, an additional important objective has been to find 
numerical techniques which can be extended to the two-dimensional 
problem. 

The paper is divided into three parts. Section I treats the physics of 
the model and Section III gives numerical results. Section II deals 
briefly with the numerical techniques used and it may be omitted by 
the reader without loss of continuity. 
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Fig. 1 — Buried-channel device structure. 



The buried-channel device has a layered structure wjiich will be 
modeled in one dimension as follows (Fig. 1): ' . 

£ * ^ hi', silicon dioxide (Si0 2 ) with relative dielectric constant 
ei/e =s 4. (e D is capacitivity of free space.) 

hi ^ x ^ h 2 : p-type silicon with acceptor number density Wa(x) 
and relative dielectric constant e 2 /t = 12. : 

h 2 ^ x ^ oo : n-type silicon uniformly doped with constant donor 
number density Nt> and dielectric constant e 3 — *i. 

The point x = is a perfectly conducting boundary held at a potential 
V . The potentials 4>i(x), faix), and 4>z(x), in the-Si0 2 , p-type\ and 
n-type regions, respectively, satisfy the one-dimensional Poisfcon 
equation, 



d*<t>i{x) 



= -p.oo/e.v t = 1,2,3, * ;•; \<i): 



dx 2 

where the pi(x), the lineal charge densities, are nonlinear functions of 
the potentials <tn{x). ; ,.' ■ 

The functional forms of the pi(x) are determined by' the assumptions: 

(i) The Si0 2 is a perfect insulator. 
(ii) The generation and recombinatioh rates for hole's and electrons 

are zero. > • . 

(in) Hole and electron currents are zero at the time of observation. 
(iv) The flat-band voltage is zero. 

(v) The injected free holes and electrons are separately in 
equilibrium. ' . j 

Conditions (i) through (v) define a static device for which the p.-(af) are: 

Pi(z) =0, "I ■ . i 

P2(x) =qlp(x)-N A (xn <\: ■ (2) 

Pi(*) = <£#» + *»(*) T *MlJ . 

where w(x') and pCx) are the number densities df free electron^ and 

holes, respectively. The most general expression' for rtt» 3B 

n(z) = /"" gimrXfydti,^ ' .' , • (3) 
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where g{E) is the density of states, F e (E) is the Fermi distribution for 
electrons, and E c {x) is the conduction band edge. Substitution for g(E) 
and F.(E) yields 

n{X) * V .. <., 1 + exp «E - E,)/kT) (4) 

This specific choice for the density of states corresponds to the simplest 
possible band structure. E/ is the equilibrium Fermi level for electrons 
and N° c is the constant 

N° c = 47r(2m c /h 2 )*, (5) 

where m c is the effective mass of an electron in silicon and h is Planck's 
constant. 6 

A similar expression for the distribution of holes is 

P(x) =( M Z) g(E)F h {E)dE, (6) 

J — CO 

where F h (E) = 1 - F,(E) ; substitution for g(E) and F h (e) yields 

p(x) = N° f ™ <*■<*> ~ E)idE • (7) 

m) ff 'L 1 + exp «E fh - E)/kT) (7) 

E v (x) is the valence band edge and E fh is defined as the pseudo-Fermi 
level 7 for the holes. N° is the constant given by 

N° = 47r(2m„/h 2 )*, (8) 

where m v is the effective hole mass in silicon. 8 E c (x) and E v (x) are func- 
tions of <t>i(x) given by: 

E e (x) = E°- qfrix), (9) 

E v (x) = E% - q<t>i(x). (10) 

E° and E° c are the valence and conduction band edges at x = « , and 
q is the magnitude of an electronic charge. p(x) and n(x) are functions 
of the yet undetermined constants E f and Efh. Later they will appear 
only in the forms 

v = (E°- E,)/kT, 
7,' = (E° - E fh )/kT. 

Since only difference terms appear, no energy reference need be estab- 
lished for the model, y will be determined from the electron charge 
neutrality condition at x ™ oo . q' is fixed by the total amount of free 
positive charge Q+ = f hl °° p(x)dx in the device. In reality, determina- 
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tion of r\' is difficult, so the alternative scheme of choosing v' and cal- 
culating the resulting Q+ is employed. Substituting eqs. (2), (4), and (7) 
into eq. (1) results in the system of differential equations which charac- 
terize the device described above. In the appendix these equations are 
reduced in a straightforward manner by using the Boltzman approxi- 
mation; justification for its use is given. The resulting equations which 
will be solved are (' denotes d/dy) : 

flXy) = 0, ^ y ^ h/\, (11) 

tf(y) = a - dm,/mJ*e"r+*M, hi/\ ^ y ^ h*/\, (12) 

t!(y) - e* tiv) - 1 - (m v /m e )*e'"-+ tW , h*/\ £ y ^ », (13) 

where 

4, = q <j>/kT, 

y = x/\ 

x = VkZVgWjy, 

N t = N°(kT)*, 
N c = N° c (kT)i, 

t a = the bandgap energy in kT's, 
v = (E° — E f )/kT is a constant dependent 

on doping levels, 
v ' = (E° c — E fh )/kT is a parameter depend- 
ing on the pseudo-Fermi level, 

p = V + v' — e a , 
a = N A /N D . 

m, and m e are the effective masses of holes and electrons in silicon; 
they are 1.08m and 0.59 m , respectively. 9 Equations (11), (12), and 
(13) satisfy boundary conditions 

*i(0) = V (14) 

and, consistent with eqs. (9) and (10), 

^ 3 (°°)=0. (15) 

At y = hi/\ and y = hi/\ the continuity conditions are 

*i(Ai/X) = Hhi/X) (16) 

Hh2/\) = ^(VM (17) 
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II. NUMERICAL SOLUTION 

The system of differential equations (11), (12), and (13) together 
with boundary equations (14), (15), (16), and (17) are solved numeri- 
cally by the method of finite differences. 

For computational convenience the length of the device will be 
truncated from the whole half-line to the segment [0, L]; the distance 
L is chosen such that | ih( °° ) — ^i{L) | and | ^ 3 '( <*> ) — &a(L) \ are 
suitably small. [0, L] is partitioned by a mesh of N points and the 
solution at each point y, is denoted by \f/\ The mesh lengths (distance 
between successive points) are 5 , 8 P , and S„ in the oxide, p-region, and 
n-region, respectively, and are chosen such that Vni = hi/X and 
y n, = h 2 /\. The boundaries are ^° = ^i(O) = V„ and ip N = f 3 (L) 
« ^ 3 (oo) =0. The subscripts on \f/(y) can now be dropped since the 
superscripts identify the solution point. The second derivative is ap- 
proximated by the second difference 

<W?/) _ ,,,,, O.,,- , ,HW« 



dy'< 



= (+ i+l - 2f< + t*- l )/S* 



in each region; 5 is one of 8 OI 8„, or 8 P as appropriate. Using this ap- 
proximation to discretize eqs. (11), (12), and (13) results in the matrix 
equation 

MO<]= IXy«,*9] (18) 

which has rows generated by eq. (19); p(yt, $*) is defined in the three 
regions by the right-hand side of eq. (19). 

f0, i < N h 

p+i - ty + ^•- 1 = U<r - (m,/m e )*e'*-* t )8 v , N x <i< N a , (19) 
[( e W _ 1 _ (m,/w» e )*e"«-**)fi BJ i > N 2 . 

Rows corresponding to the solutions at the boundary points t/at, and 
yti, are obtained from eqs. (16) and (17) where the first derivative is 
approximated by 

=Ff(y) « (+W(y) - W(V±8) 

+ fy(y ± 25) - 2+(y ± 35))/65, (20) 

and 8 is the appropriate mesh size for each region (the positive or for- 
ward derivative uses points to the left of y). Equation (20) results from 
simultaneous solution of the Taylor expansions for \f/(y ± 8), \f/(y ±25), 
and \J/(y ± 35) with terms 0(5 4 ) and higher dropped; the result is ac- 
curate to 0(5 3 ) which is consistent with the second difference approxi- 
mation to d 2 \J//dy 2 . In the oxide layer, a first-difference approximation 

fix) = (*» - **-*)/«. (21) 
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is adequate since the solution in this region is linear. The boundary- 
equations are: 

for i = Ni, 

-(6* p ci/*.€»)^ w »- 1 + (11 + 6W«o€2)* Wl 

- W" + fy Nl+i - 2^+ 3 = 0; (22) 
for i = N 2 , 

(-1/*p)(2^ jv »- 8 - 9^'- 2 + 18^- J ) + ll**»(l/*» + 1/« P ) 

- (1/S n )(18^ s+1 - fy N > +2 + 2^ Ar »+ 3 ) = 0. (23) 

The matrix A in eq. (18) is tridiagonal except for the Nith. and JV 2 th rows. 
For each choice of the parameter p , the [$*] of eq. (18) are solved for 
iteratively using a combination of successive under- and over-relaxa- 
tion 10 (SOR) on the equation 

[>•'] = M-»[>&**0J ( 24 ) 

The procedure described here differs slightly from the usual SOR in 
that the transformed equation (24) is solved instead of eq. (18). For 
the t'th row of eq. (24), the j -f- 1th estimate of ^* is given as 

fti = # + "(-A*' - #), (25) 

where w is a relaxation parameter with values < o> ^ 2. $' is a New- 
ton's method solution of eq. (26), a transcendental equation in ^result- 
ing from the ith row of eq. (24) with the remaining N — 2 variables ^* 
held constant. The coefficients Out 1 in eq. (26) are elements of A' 1 . 

f - E a* l p{y k , ftO - a^piyi, f ) - "£ a£ l p(Vk, *5) = 0. (26) 

A-l k-»+l 

Criteria for convergence of the process, as well as the choice of the 
value of w, is based on the residual rj defined at the jth. iteration as 

r,= Z |*}-fti I • (27) 

It can be shown that if (18) were a linear system of equations and if f 
were the true solution at ?/, then 

sup If -ft.il/aq> \?\*CMr M (28) 

where C(co) is a constant dependent on the choice of »." For the optir 
mum value of v, C(oj) is 0(N) while for values of cu only slightly different 
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from the optimum C(oj) can be 0(iV 2 ). In the computations presented 
in the next section, N ^ 100 so the iterative process was stopped when 
the residual terms were less than 10 -8 ; this allowed margin for the fact 
that (18) is nonlinear. A discussion of the choice of to is necessarily even 
more heuristic. It was found that for a mesh of N = 50 points and large 
negative values of p corresponding to small positive charge densities 
(i.e., e"»~*'' « 1 for all i), the iterative scheme was convergent for any 
choice of cu and any initial estimate of the {^*} . For less-negative values 
of p and the corresponding larger values of e Po "~*', the scheme was con- 
vergent for over-relaxation (w > 1) only if r, < 10 for all j (approxi- 
mate figure) ; but using under-relaxation the process was well behaved 
with r/s as great as 10 4 . It should be pointed out that due to the ex- 
ponential nature of the right-hand side of eq. (19) even a very good initial 
estimate typically resulted in ri > 10 3 when over-relaxation was ap- 
plied. Although the SOR process described above may well be stable 12 
regardless of r, values, the magnitude of the exponential terms in (19) 
limits machine computations. The combination of under- and over- 
relaxation detailed below eliminates this problem. 

Good initial estimates for the {^'} were obtained by choosing the p„ 
Values with equal spacing Ap; the first p value being the small positive 
charge-density case described above. With Ap < 10 a linear estimate 
of the {$'] Tor successive values of p was adequate. The iterations were 
under-relaxed (« = 0.5) until the residual term was less than 10; suc- 
cessive iterations were over-relaxed so as to increase the rate of con- 
vergence. The choice of as for the SOR steps was not critical; values in 
the range 1 < w 2g 1.5 had approximately the same rate of conver- 
gence. Values above 1.5 did not converge for all values of p . 

Total positive charge in the device for a given value of p is calcu- 
lated by 

Q + = f ""' e'°-**Mdy + r 6 e"-**^dy. (29) 

Using a piecewise linear approximation for \J/(y) on [ijk, Vk+\], 

+(y) SS P + ((» - ¥*)/*)(#*» - **), (30) 

and summing the integrals over each such interval in [_Vn v L~] gives the 
approximation 

Q + ^S P £ e-ie-*' - e-* i+l )/(f i+1 - *»') 

+ S n Z e*>(e-*< - e-*' + W' +1 - M, (31) 
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which can be modified to include the case \f/ i+1 = ip\ Dimensional 
charge in coulombs can be calculated from Q+ by Q = Q+Nd^- 
The operational scheme may be summarized as follows: 

(i) Choose a value of p corresponding to small total charge Q + 

( Po -» - * in eq. (24) => Q+ -^ 0), 
(it) Solve for \f/(y), thus determining p(y), 
(Hi) Integrate p(y) to find Q+, 
(iv) Increment p by Ap and repeat (ii), (in), and (iv). 

The technique of starting with Q+ « and slowly adding positive 
charge to the device is important since it is this scheme of operation, 
in conjunction with the successive under- and over-relaxation, that 
avoids the exponential overflow limitations in machine computation. 

III. COMPUTATIONAL RESULTS 

In this section, numerical results for two specific device configura- 
tions will be given. The first device has a constant doping profile 
(Na = a constant) and dimensions 

hi/X = 0.48, 
h 2 /\ = 12.48, 
L = 42.48. 




4 8 12 16 20 24 28 32 36 40 44 

x/\- DISTANCE IN DEBYE LENGTHS 

Fig. 2 — Potential versus distance for a buried-channel device with uniformly 
doped p-region and V„ = — 4 volts. 
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Fig. 3 — Charge density versus distance for a buried-channel device with uniformly 
doped p-region. 



X in this case is approximately 0.415 /mi. Doping levels are 

N A = 2 X 10 15 cm" 3 , 
JV D = 1X 10 14 cm" 3 , 
so 

a = Na/Nd = 20. 
N e is calculated by 13 

N c = 4.831 X 10"(m«/m.)»r», 

with T = 300° K and m c as before. 
Using eq. (46), 77 is found to be 

77 = 12.09. 

The nondegeneracy condition for this device may now be stated using 
eq. (45) 

6, - 17' + +(y) > 3.5, 
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Fig. 4 — Potential at the Si-Si0 2 interface versus total positive charge in a buried- 
channel device with a uniform p-region and V = — 4 volts. 



or 



V' ~ e g - +(V) < -3-5. 

Adding rj to both sides of the inequality gives 
Po - Hy) < V - 3.5. 




-800 
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Fig. 5 — Potential versus distance for a buried-channel device with a Gaussian 
p-region doping profile and V = —4 volts. 
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Fig. 6 — Charge density versus distance for a buried-channel device with a Gaussian 
doping profile. 



Then for each choice of the parameter p , the solution set {^'} must 
satisfy: 

Po - +< < 8.59. (32) 

From eq. (12) it is clear that if one assumes the Boltzmann approxi- 
mation to hold, then as long as Q + < <r(h 2 — hi)/X the equilibrium con- 
dition causes the lineal charge density to have constant sign so 



Po — f* < log <j(m c /m v )* 



(33) 



for all i; for this device, eqs. (32) and (33) are always consistent if 
N A < 1.3 X 10 18 cm- 3 . 

Computation was performed using a 90-point mesh and took ap- 
proximately 4.8 minutes of processor time on a Honeywell 6000-series 
machine. 

Figures 2, 3, and 4 summarize the computational results. Figure 2 
shows potential solutions {if/*} versus the points ?/, for a family of p„ 
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Fig. 7— Potential at the Si-Si0 2 interface versus total positive charge in a buried- 
channel device with a Gaussian p-region doping profile and V a = — 4 volts. 

values with V = — 4 volts. The range of total charge values is indi- 
cated. Figure 3 is a plot of the linear charge density with Q+ values 
indicated. Figure 4 is Q+ plotted against $ Nl (the oxide-interface 
potential). 

A realistic modification to the device studied thus far is to allow a to 

have y variation. The second set of results presented are for the same 

device as described above but with the p-region having a doping profile 

a(y) = (2>. + l)tf-l<r-»i>/(ta-ii)l« ^ (*>■+» _ 2. (34) 

The values for h h h 2 , and V are as before, L = 32.48, and D, = 46. 
[This value of D, corresponds to an average doping level in the p-layer 
of N A = 1.9 X 10 16 /cm 8 , and eq. (34) describes a doping profile as if 
the p-layer were formed by drive-in diffusion. 14 ] The solutions \\l/ 1 } 
must still satisfy eq. (32). Figures 5, 6, and 7 are a summary of results. 
Comparison of the two cases shows that for the doping profile in eq. 
(34) the potential minimum is greater and the "channel" is shifted 
toward the oxide. In both cases the added positive charge is contained 
entirely in the p-region and it fills the region starting from the side 
remote from the oxide interface (see Figs. 2 and 5). 
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APPENDIX 

Substituting eqs. (2), (4), and (7) into eq. (1) results in (' denotes 
d/dx) : 

4>'i(x) = 0, 0gx£h h (35) 

Mx)--]N A -N.J_ m 1 + exp L{Efh _ E)/kT1 j- , 

h ^ x ^ hi, (36) 
»/ x g|»„r (E-E e (x))*dE 

m f '•« (E 9 (x) - E)»dE \ 

"•}-* 1 + exp [(£,„-£) AT] J' 

ht ^ x S w . (37) 

Transform the integral involving E c (x) by making the substitutions: 

e = (E - E c (x))/kT, 
E c (x) = E° c - q<p{x), 
v = (El - E f )/kT. 

For the integrals involving E v (x) make the substitutions: 

e = (E v (x) - E)/kT, 

E,(x) = E% - q V (x), 

rj' = (E° c - E fh )/kT, 

e = (E° - E°)/kT. 
For both cases let 

4>{x) = q<p(x)/kT 

and make the change of variables 

V = */\ 
where 

x = VkTVgWD. 

These substitutions and some straightforward manipulation reduce 
eqs. (35), (36), and (37) to 

tib) =0, ^ y g hi/\, (38) 

,»(»\ __N A _N% n-n,/"* e^de 

* 2{y) - N D N D {kI) in 1 + exp (e + «. - r,' + fcfo)) ' 

fej/X ^ y ^ A 2 /X, (39) 



1022 THE BELL SYSTEM TECHNICAL JOURNAL, JULY-AUGUST 1973 

tfdr) = jt d (kT)*J o 1 + exp p + - _ ^ |(y) -, - 1 

tf D ^ J 7o 1 + exp (£ + «.- f,' + **(»)) ' 

A2A ^ 1/ ^ ». (40) 

These equations are further simplified by making the classical Boltz- 
mann approximation for the electrons : 

r *** ~r ^e (41) 

Jo 1 + exp (e + v ~ t(.y)) Jo exp(e + 77.- $(y)) 

which holds since y - \p(y) is always positive by several kT (see Fig. 
8). The last result is integrable and simplifies the approximation to 

T(f) exp (+(y) ~ v), (42) 

where r(-) is the usual gamma function. This approximation has been 
shown to be less than one percent in error if 16 

t, - +(y) > 3.5. (43) 

Similarly, for the holes, 



/ 



€^€ 



^ r(|) exp w - i„ - +(y)) (44) 



o 1 + exp (e + e - y' + *(y)) 

if 

e B ~ v' + *(V) > 3.5. (45) 

The validity of this last approximation is not clear since +(y) is de- 
pendent on i}' and \p(y) becomes large and negative. Note that requiring 
the inequalities (43) and (45) to hold is equivalent to requiring the 
device to be nondegenerate. Computational results show this approxi- 
mation to be consistent. 

77 is expressed in terms of the constants N c , N v , and N D by requiring 
electron charge neutrality at x = -» . It should be stated that the cor- 
rect condition to use at this point is complete charge neutrality. This 
requirement, however, is computationally indistinguishable from the 
algebraically simpler condition used here. Using eqs. (40) and (15), the 
neutrality condition is 

(N e /N D )e-"T(%) -1=0, 

so 

N C /N D = e'/m). (46) 
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Fig. 8 — Band diagram for buried-channel device. 

Equation (46) defines the Fermi level for the electrons. Using eqs. (5) 
and (8), 

N V /N D = (N,/N C )-(N C /N D ) = (m./w.)Wr(t), (47) 

where m v and m c are the effective masses of holes and electrons in 
silicon; they are 1.08 m and 0.59 m , respectively. 10 Using eqs. (42), 
(44), (46), and (47) in eqs. (38), (39), and (40) results in equations 
(11), (12), and (13) which are to be solved. 
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